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We describe a criterion for particles suspended in a randomly moving fluid to aggregate. 
Aggregation occurs when the expectation value of a random variable is negative. This 
variable evolves under a stochastic differential equation. We analyse this equation in 
detail in the limit where the correlation time of the velocity field of the fluid is very 
short, such that the stochastic differential equation is a Langevin equation. 



1. Introduction 

1.1 Illustration and context 

Figure ^ illustrates a simulation of the distribution of small particles suspended in a 
three-dimensional random flow: the particles are modelled as points, but they are shown 
as small spheres to make them visible in the figure. The suspended particles do not 
interact (that is, the motion of each particle is independent of the coordinates of the other 
particles, the equations of motion arc given below). We show the initial configuration 
(Figure ^i), and two snapshots of the particle positions after a long time, with differing 
values of the fluid viscosity, (Figures^ 3 an d c). fn one case the particles aggregate, in the 
sense that the trajectories of different particles coalesce, fn the other their distribution 
shows some degree of clustering, but their trajectories never coalesce. In this paper we 
present an analysis of the transition between aggregating and non-aggregating phases, 
which we term the 'path-coalescence transition'. 

There are numerous experimental observations that when small particles are suspended 
in a complex and apparently random flow, their density becomes non-uniform. Clustering 
of particles into regions of high density has been observed in e xperiments on particles 

floa ting on the sur face of liquids with a complex or turbulent flow ( Sommcrer fc Ott 199.lt 

ICressman et al. 20041) . and also turbulent flow in channels l|Fessler et al 1994tlvan Haarlem et al. 1 9981. 

The conditions under which this occurs are not yet fully understood. The aggregation 

effect illustrated in figure 1 is an extreme form of clustering. There appears to be less 

experimental work on aggregation, but coalescence of su spended wat er droplets is clearly 

very important in the formation of rain drops in clouds (jShaw 20031) . Even less is known 

about aggregation. In particular it is not clear when a model of particles suspended in 

a random flow can exhibit aggregation, and when additional physical phenomena (such 

as differential drift velocities under gravitational forces, or Brownian diffusion) must be 

invoked. 

Earlier theoretical work on clustering in random flows has used Fokker-Planck equa- 
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Figure 1. Illustrating the aggregation of non- interacting particles in a random 
three-dimensional flow: the motion is defined by equations 1)1. 1^ to l|1.3jl . (using a simplified 
form explained in section [3] see equation 13.231 1. The left panel shows the initial configuration 
at t — 0, the center and right panels show final configurations (at t = 180 r) for two different 
values of 7. In the case shown in the center panel, trajectories coalesce, until eventually all 
particles follow the same trajectory. The particles in the right panel exhibit density fluctuations, 
but the trajectories do not coalesce. 



tions determining moments of the particle density of advected particles l|Klvatskin fc Gurarie 199Sh . 
This 'passive scalar' approach does not allow for the effects of inertia of the particles. It 
has been supplemented by a pertur bative analysis of the deviations of particles from ad- 
vected traject ories, as proposed b v lMaxev (1987ll. These two a pproaches are combined 
in papers by lElperin et aJ (199fj) and iFalkovitch. et. al (200111 . Numerical work indi- 
cates that clustering in solenoidal flows occurs most readily when the correlation time of 
the velocity field is comparabl e with the time constant associated with the viscous drag 
ijSigurgeirsson fc Stuart 2002j) . The aggregation (as o pposed to cluste ring) of particles 
by random flows has received relatively little attention. iDeutsch (1985|) appears to have 
been the first to propose that particles subjected to a smooth random flow can coalesce, 
and showed numerical evidence that this can happen in one dimension. He argued that 
there is a transition between coalescing and non-coalescing phases, and identified the 
dimcnsionless parameter which determines the phase transition in one dimension. 

This paper describes results obtained from a new approach to characterising particle 
aggregation in random flows. It is based upon calculating a Lyapunov exponent describ- 
ing the rate of separation of nearby particles from a solution of a system of stochastic 
differential equations: the Lyapunov exponent is the expectation value of one of the 
variables. In the limit as the correlation time of the flow approaches zero, the stochas- 
tic differential equations can be reduced to a pair of Langevin equations. Our results 
for three-dimensional flows which are di scussed here build upon ea rlier work by two 
of the present authors in one dimension (Wilkinson & Mehlig 2003) (whe re we solved 
Deutsch's model exactly) and two dimensions l|Mehlig fc Wilkinson 20041) . The three- 
dimensional case is most important for physical applications, but involves substantial 
additional technic al complications. 

We remark that lPiterbarg (200l[) also considered the two-dimensional case, and quotes 
an analytic expression for the maximal Lyapunov exponent. His expression is incorrect 
in two dimensions, because the generating function that he uses is divergent for the 
equilibrium distribution. His calculatio n can be adapted to give the correct expression in 
the one dimensional quoted in lMehlig fc Wilkinson (200411. 

The definition of Lyapunov exponents is explained bv lEckmann fc Ruelle f 1979() , who 
also discuss the method we use for extracting the Lyapunov exponents from direct numer- 
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ical simulations. We note that in the case where inertia of the particles can be neglected, 
our results reduce to calculating the Lyapunov exponent for a spatially c orrelated Brow- 
nia n motion, which was discusse d from a mathematical point of view bv lLe Jan (19851) 
and lBaxendale &: Harris' 



1.2 The model 

The most natural model for theoretical investigation is the motion of spherical particles 
(radius a, mass to) moving in a random velocity field u(r, t) with specified isotropic and 
homogeneous statistical properties. The particles are assumed not to affect either the 
flow, or each other's motion, and to experience a drag force given by Stokes's law: the 
force fdr on a particle moving with velocity v relative to the fluid is fdr = —Qirrjav, where 
•q is the viscosity of the fluid. We will simplify the problem by assuming that the particles 
are made of a material which is much denser than the fluid in which they are suspended: 
this enables us to neglect the inertia of the displaced fluid. Accordingly, we consider a 
large number of suspended particles, initially with random positions and zero velocity, 
having equations of motion 

r=— p, p = -7[p- mu(r,t)] . (1.1) 

TO 

The random velocity field u(r, t) could be either externally imposed (for example, if a 
gas is driven by an ultrasonic noise source), or self-generated (as in the case of turbu- 
lence). The equations of motion with displaced mass effects included are discussed by 
lLandau fc Lifshitz (1958). Our neglect of displaced mass effects is justified for aerosol 
systems. 

In order to fully specify the problem we must define the statistical properties of the 
random velocity field u(r,t). The random force f = TO711 is generated from a vector 
potential A = (A±, A2, A3) and a scalar potential <f> = Aq: 

f = V<£ + VAA. (1.2) 

The scalar fields Aj(r, t) have isotropic, homogeneous and stationary statistics. We as- 
sume that these fields are statistically independent, and that they all have the same 
correlation function, except that the intensity of <f> — Ao exceeds that of the other fields 
by a factor I /a 2 : 

(A i (r,t)A j (r',t'))=S ij {(l-a 2 )S i0 +a 2 ]C(\r-r'\,t-t') . (1.3) 

The random force f(r, t) is characterised by its typical magnitude a, and by the corre- 
lation length £ and correlation time r of the correlation function C(R,t). In the case of 
a well-developed t urbulent flow , the velocity field has a power-law spectrum with upper 
and lower cutoffs IjFrisch 1997(l . If a random velocity field modeling fully-developed tur- 
bulence is used in our model, it is most appropriate to take the correlation length and 
time to be those of the 'dissipation scale', that is, the cutoff with the smaller length scale. 

The system of equations is characterised by three independent dimensionless parame- 
ters, which we take as 

v = IT , x = — 7 , a ■ (1-4) 

to£ 

The parameter v is a dimensionless measure of the degree of damping, x is a dimensionless 
measure of the strength of the forcing term and a measures the relative magnitudes of 
potential and solenoidal components of the velocity field (which is purely potential when 
a = 0, and purely solenoidal in the limit as a — » 00). 
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1.3 Description of our results 

We show that the phase transition is determined by a Lyapunov exponent, Ai, de- 
scribing the separation of nearby particles: their trajectories coalesce if Ai < 0. Here we 
describe a new and general approach to this problem, reducing the determination of the 
Lyapunov exponent to the analysis of a simple dynamical system, described by a sys- 
tem of ordinary differential equations containing stochastic forcing terms: the Lyapunov 
exponent is found to be proportional to the expectation value of one coordinate. These 
stochastic differential equations are derived in section They introduce an apparent 
paradox: the structure of the equations appears to be identical in two-dimensional and 
three-dimensional flows, suggesting that the path-coalescence transition is fundamentally 
the same in two and three dimensions. That would be a very surprising conclusion. 

In order to demonstrate and illuminate our method, in sections [21 and 0] we pursue 
the solution of this problem in considerable depth in one limiting case, namely the limit 
where the correlation time r of the random velocity field is small and the random force 
is sufficiently weak (strictly, we consider the case where \ "C v <C 1). In this limit, 
the system of ordinary differential equations described in section [21 becomes a system 
of two coupled Langevin equations. In section [21 we show that there is in fact a dif- 
fer ence between the three-dim ensional problem and the two-dimensional case studied 
in iMehlig fc Wilkinson (2 004) , which is a rather su btle example of the difficulties in 
applying Langevin approaches to nonlinear equations Ijvan Karnnen 1 992). 

In section 0] we discuss a perturbation theory for the Lyapunov exponent describing 
the phase transition, expanded in powers of a parameter e = \ v ~ 3 ^ 2 which is a dimen- 
sionless measure of the inertia of the particles. The perturbation theory is constructed by 
transforming the Langevin equation first into a Fokker-Planck equation, and then into a 
non-Hcrmite an perturbation of a three-dimensional isotropic quantum harmonic oscilla- 
tor. We are t hen able to use the harmonic-oscillator creation and annihilation operators 



( Dirac 19301) to express the perturbation theory in a purely algebraic form, enabling us 



to compute the coefficients to any desired order. We investigate the phase diagram (the 
line in parameter space separating coalescing and non-coalescing phases), using both 
Monte Carlo averaging of the Langevin equation and the results of our perturbation 
theory. We find that aggregation can only occur if the random flow has a certain de- 
gree of compressibility, which increases as the effects of inertia increase, until there is no 
coalescing phase even for a purely potential flow field. 

The analysis in sections [21 to 0] considers the case where all of the suspended particles 
have the same mass m and damping rate 7. This ideal can be approached quite accurately 
in model experiments, but in most applications in the natural world and in technology, 
the suspended particles will have different masses, sizes and shapes. In section [S] we 
discuss the effect of dispersion in the distribution of masses in the one-dimensional case: 
these arguments can be adapted to treat higher dimensions and dispersion of the damping 
constant 7. We argue that path coalescence is not destroyed by mass dispersion (although 
of course it is no longer a sharp transition) . In this paper we give a quite comprehensive 
discussion of the the case where the correlation time of the random flow is very short, 
and the stochastic differential equations derived in section 2 can be approximated by 
Langevin equations. Section 6 discusses how our approach can be extended to other 
cases. 
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2. Equations determining the Lyapunov exponent 

To determine whether particles cluster together, we consider two nearby trajectories 
with spatial separation 5r and momenta differing by <5p. The linearised equations of 
motion derived from are 

Sr= — Sp, Sp = -jSp + F(t)Sr . (2.1) 
m 

Here F(t) is proportional to the strain-rate matrix of the velocity field, with elements 

F l3 {t) = ^±(r(t),t)=rn^(r(t),t) . (2.2) 

It is convenient to parameterise <5r and 5p as follows: 

Sr = Xm , Sp = X(Ymi + r 2 n 2 ) , (2.3) 

where ni and n 2 are orthogonal unit vectors, which depend upon time. The parameter 
X is a scale factor: trajectories coalesce if X decreases with probability unity in the long- 
time limit. In the three-dimensional case, we find it convenient to introduce the third 
element 113 = ni A n 2 of a time-dependent orthonormal basis so that 

rii.rij = Sij and = e^n,- A n fe . (2.4) 

Differentiating (|2.3(l . and substituting the resulting expressions into H2.lt gives 

Sr = ini + Xhi 

= -X(Y 1 n 1 +Y 2 n 2 ) (2.5) 
m 

dp = XiY.m + Y 2 n 2 ) + X(Fmi + Y 2 n 2 ) + X(y ini + Y 2 h 2 ) 
= -7X(y ini + y 2 n 2 ) + F(t)ni . 

Projecting Sr onto the unit vectors gives the following three scalar equations of motion 

ni.Sr = X = —Y X X 

m 

n 2 .6r = X(n2.rn) = -Y 2 X (2.6) 
m 

n^.Sr = Xhi.ris = . 

The last of these equations implies that rii A n 2 = 0, so that ni is proportional to n 2 : 
we write 

ni = 9n 2 (2.7) 

so that the equation for n 2 .5r gives 

e=-Y 2 . (2.8) 
m 

The first equation of l|2.6|l indicates that X is a product of random variables, and therefore 
has a log-normal distribution, that is, the logarithm of X has a Gaussian probability 
density. In the limit as t — > 00, the mean and variance of log c X are both linear functions 
of time: 

(log e X) - X x t + ci , var(X) = fit + c 2 (2.9) 

where Ai, fx, c\ and c 2 are constants. If Ai < 0, the probability of log X exceeding any 
specified value approaches zero as t — > 00, implying that trajectories of nearby particles 
almost always coalesce. 
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The Lyapunov exponent Ai is the mean value of the derivative d log e X/ dt, so that the 
first equation of (|2.6J) gives 

Ai = -(Yi> . (2.10) 
m 

Now consider the three projections of <5p, as given by equation (j2.5(l : 

ni .(5p = XY t + XY t + XY 2 (n 1 .h 2 ) 

= -~iXY x +m.F(f)niX 
n 2 .£p = XY 2 + XY 2 + Jfyi(ni.na) 

= -jXY 2 + n 2 .F{t)mX (2.11) 
n 3 .<5p = XYi(n 3 .ni) + XY 2 (n 3 .h 2 ) 

= n 3 .F(t) ni X . 

We introduce the notation 

n l {t).F(t)n. ] (t) = F' l3 {t) (2.12) 

and note that the statistics of the transformed matrix elements F-j(t) are the same 
as those of the original elements Fij(t), because the statistics of the velocity field are 
isotropic. Using eqs. I|2.6|l to l|2.8|l and (rij.nj) + (ri,.iXj) = to simplify, we find the 
following equations of motion for the variables Yi 

Yi = -jYi + -(Yi - Y?) + Fi^t) (2.13) 
m 

Y 2 = - 1 Y 2 --Y 1 Y 2 +F! 21 {t) . 
m 

Finally, the equation for n 3 .(5p gives 

n 2 .n 3 = ■ ( 2 - 14 ) 

Eqs. 12.10|l and l|2.13|l are the principal results of this paper. Eq. I|2.1U|I shows that the 
Lyapunov exponent (the sign of which determines whether or not path coalescence occurs) 
is given by the expectation value of a random variable Y\ of a simple, finite dimensional 
stochastic dynamical system, described by eqs. (I2.13|) . This dynamical system is almost 
completely de-coupled from the other variables: the equations for Y\ and Y 2 do not 
depend upon X, and the vectors rij(i) only enter these equations through the evaluation 
of the random matrix elements F^ . We pointed out that statistics of these elements are 
independent of the orientation of the orthogonal triplet (ni, n 2 , n 3 ). 

In the two-dimensional case the analy sis leading t o (12.1311 proceed s along similar lines, 
and leads to the same pair of equations l|Mehlig fc Wilkin son 20041) . The only difference 
is that the equation (|2.14|l is absent in the two-dimensional case. This suggests that the 
expression for the Lyapunov exponent Ai = (Yi) jm should be the same in two and three 
dimensions. This would be a surprising conclusion, but it is not obvious how it can be 
averted. However, it does prove to be false, as will be demonstrated in the next section 
for the limiting case where xC^Cl- 



3. The Langevin approximation 

Let us now consider how to treat eqs. I|2.13|l in the limit where the correlation time r of 
Fit) is very short. Because the random field f (r, t) is fluctuating very rapidly, the position 
r(t) of a particle at time t is independent of the instantaneous value of the force f(r, t), 
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so that the value if Fij(t) = dfi/dxj(r(t),t) at the position of the particle is statistically 
indistinguishable from a random sample of the field dfi/dxj. We also assume that the 
gradients of the fluctuating forces (the quantities F[ x and F 22 in equations H2.13fl ) are 
sufficiently small that the typical magnitude of the displacement of the variables Y%, Y 2 
occurring during the correlation time r is small compared to the typical magnitude of 
these variables. This condition is expressed in terms of the dimensionless variables \ 
and v later in this section. Under these conditions of short correlation time and small 
amplitude, equations (|2.13|) may be replaced by Langevin equations. At first sight, this 
would appear to lead to 



dYi = 
dY 2 = 



_ 7 y 1 + l(y2_ y2) 
m 



- 7 y 2 - —y x y 2 

m 



dt + d/i 
dt + d/ 2 



(3.1) 



where the dfi are increments of a Brownian process, satisfying (d/i) = and (d/jd/j) = 
2Dijdt, for some const ant diffusion coefficients D g. In the two-dimensional case, this 
expectation is correct l|Mehlig fc Wilkinson 2004j) . and eqs. ijSJJ are the appropriate 
Langevin equations. In the three-dimensional case, we will see that an additional drift 
term must be added to the second of these equations. This is a consequence of the fact 
that, in the three-dimensional case, Y 2 is a non-linear function of the components of the 
vector (5p — Y\b~v because it is the magnitude of this vector. This is an example of the 
difficulties t hat arise when treati ng Langevin equations involving non-linear functions of 
noise terms ijvan Kampen 1992]) . 

In order to determine the correct Langevin equations to model (|2.13|l . let us consider 
the integral of the stochastic forcing terms dfi over a time interval 8t which is long 
compared to t, but short enough that the change in the variables Yi occurring in time 
5t can be neglected. We define 



and find 
where 



/t+st 
dt'i^(i') 



{8f i 5f j ) = 2D ij 5t + 0{r) 



D = i 

U%] — 2 



dt (0)) . 

The calculation of (Sfi(t)) is more subtle. We have 

(Sfi) - / dt (n i (t).F(t)n 1 (t)) 



We must take account of the fact that the unit vectors iii(t) are rotating: we write 

3 

fc=i 

where Rik (t) are elements of a rotation matrix. We now write l|3.5|l in the form 

r st 3 3 



(3.2) 

(3.3) 
(3.4) 

(3.5) 
(3.6) 

(3.7) 



fe=l 2=1 
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Now consider the rotation of the unit vectors: using 1|2.10|1 and (|2.14|1 we have 
m(t) = n^O) + 6n 2 {0)t + 0(t 2 ) , 



n 2 (i) = n 2 (O)-0n 1 (O)i 



Y 



2 JO 



dt' F^(t')n 3 (0) + 0(t 2 ) , 



(3.8) 



n 3 (t) - n 3 (0) - i / df F^(t')n 2 (0) , +0(£ 2 ) 

Y 2 Jo 



so that 



We obtain 



R{t) 




6t 
1 




1 



0(t 2 ) 



(3.9) 



k 

t 



(3.10) 



This yields 



(Sfi) 
(Sh) 



dt (F{ 1 (t)+etF^ 1 (t))=0 



5i 



dt (-9F{ 1 (t)+F^ 1 (t) + 



Y 



2 JO 



dt' F'^F'^t')} 



i ["dt f dt'iF^F^t')) 

r 2 Jo JO 



Y 



D 31 St . 



(3.11) 



The Langevin equations therefore contain an additional drift term due to the fact that 
(^/2) is non-zero: the correct Langevin equations in three dimensions are 



with 



dYi = 
dF 2 = 

(dCf 



_ 7 y 1 + l(y2_ y2) 

777 



-7^2 + -TT ^1^2 

J 2 777 



dt + dCi , 
dt + d( 2 , 







(dCidO 



2D. lJ dt . 



(3.12) 



(3.13) 



The diffusion constants Dij were defined in equation (|3.4|l . Note that in two dimensions, 
however, 1)3. If) remains valid because the term arising from the rotation of 113 is absent. 

Now consider the evaluation of the diffusion constants in terms of the statistics of the 
force f(r, t). The elements of the force-gradient matrix F are 



d 2 



p.. = 

lJ dxjdxi 



(■ilk 



dxjdxi 



(3.14) 



where euk is the "Kronecker e-symbol" describing the parity of the permutation of the 
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Then define D\ = Dn, D2 = D21 = Dai, so that 

D 1 =D \1 + — J and D 2 = D Q (- + —j. (3.16) 

We introduce a more convenient dimensionless measure of the relative importance of 
solenoidal and potential fields 

and find ^ ^ T ^ 2 in the three-dimensional case because ^ a ^ 00. It is convenient 
to re-scale the Langevin equations into dimensionless form: write 

dt' = 7 dt , Xl = J^Y, , dwi = J^dd (3.18) 



and define 



p\' 2 
TO7 3 / 2 



(3.19) 



With these changes of variables, the Langevin equations become 

dx! = [-xi + e(Txl - x\)]dt' + d Wl (3.20) 
dx 2 — [—x 2 + x^ 1 — 2exiX2]dt' + du>2 

with 

(dwi) = 0, (dwidwj) = 2<%dt' . (3.21) 

Eqs. I|3.20I3.21[I must be solved to determine the expectation value of x\ in the steady 
state. The Lyapunov exponent is then given by 

Ai = ~fe{x x ) . (3.22) 

Figure [5^ compares the Lyapunov exponent obtained from a Monte Carlo simulation 
of equations 13. 201 with a direct numerical simulation of a random flow described by 
equation The Lyapunov exponents determined from eqs. (|3.2UI) and (|3.21|) for L = 

1/3, 1 and 2 are plotted as red line s. The results are compared to numerical simulations of 
<|1-1[> , using a method described in lEckmairn fcRuelle f!97flh to determine the Lyapunov 
exponent. Because we are concerned with the limit where the correlation time t is taken 
to zero, the random flow was generated using a discrete series of uncorrelated random 
impulses, acting over a small time step dt 3> t: the impulse 

r(n+l)St 

f„(r) = / dt'f(iy,i') (3.23) 

at time nSt is taken to be of the form (|1.2fl in terms of scalar fields 4> n {^) and A„(r) 
satisfying 

(0„(r)<Mr')> - °H 2 St expflr-rf /2£ 2 )<5 W (3.24) 
and similarly for A ra (r). This implies D = 3cr 2 /(2m 2 7 3 £ 2 ). 

Now we discuss the conditions under which the Langevin equations l|3.20|l and l|3.21|l 
are a valid approximation of l|2.13|l and i|2.14[l . For this purpose it is sufficient to consider 
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the one-dimensional version of equations (|3.12|1 , namely 

Y = -7Y - — Y 2 + F(t) (3.25) 
m 

(this equation appears with a different notation in IWilkinson fc Mehlig (2 003)1. The 
Langevin equations are valid provided the changes in the value of Y over the correlation 
time r is small compared to the typical values of this quantity. This criterion can obviously 
only be satisfied if the correlation time is sufficiently short that v = jt <C 1. The criterion 
also requires the stochastic force F(t) to be sufficiently weak. The deterministic part of 
the velocity, — jY — Y 2 /m, is positive in the interval from Y — —7m to Y = 0. The 
criterion on the strength of F ~ <r/£ is that the displacement over time t should be 
small compared to the width of that interval, that is \F\t <C jm. Using the fact that 
\F\ ~ <j/£, we obtain the following criteria for the validity of the Langevin approximation: 

-«1, (3.26) 
v 

For completeness, we end this section by mentioning how equations (|3.20() and 13.2111 



differ in one and -two dimensions. The one-dimensional case was considered in (Wilkinson & Mehlig 2003^ : 
the Lyapunov exponent is given by A = (Y)/m, with Y satisfying (|3.25|l . In two dimen- 
sions, as we have alrea dy remarked, the term x^ 1 is absent from the second equation of 
EEm , and \ < T < 3 ifMehlig fc Wilkinson 200^) 1 . 



4. Perturbation theory 

We now show how to obtain an asymptotic approximation for the Lyapunov exponent 
using eqs. H3.20[l and 13.211) . These equations are equivalent to a two-dimensional Fokker- 
Planck equation ijvan Kampen (T992|) ) for a probability density P{x\, X2',t'), of the form 

d t >P = DV 2 P-V.{vP)=TP. (4.1) 

Here the diffusion constant D = 1 and the drift velocity is v = (v\, V2) with components 
v\ = —x\ +t{Tx 2 , —xf) and V2 = —X2 + X2 1 — 2ex\X2- We write J- = fio + eFx, and seek a 
steady-state solution satisfying TP — by perturbation theory in e. In order to simplify 
the application of perturbation theory, it is convenient to make a transformation so that 
the unperturbed Fokker-Planck operator J-$ is transformed into a Hermitian operator. 
Rather than proceeding to the Hermitian form directly, we first map the two-dimensional 
Fokker-Planck equation to a three-dimensional equation with a rotational symmetry (we 
seek a solution which is invariant under rotation). After making this transformation, 
we find that the corresponding Hermitian operator in three-dimensional space is the 
Schrodinger operator of an isotropic three-dimensional harmonic oscillator. The pertur- 
bation analysis can then be performed very e asily, using th e algebra of harmonic-oscillator 
raising and lowering operators, described in l|Dirac 19 30). To shorten equations, we will 
use a variant of the Dirac notation scheme: in summary, functions a, b are symbolised by 
vectors \a), |6), linear operators are denoted by a 'hat', e.g. A, and the integral over all 
space of the product of two functions is denoted by the inner product (a|6). 

In the original form, the action of the unperturbed part of the Fokker-Planck operator 
on a function P is 



F a P = (dj + d 2 )P + d 1 {x 1 P) + d 2 [(x 2 - x^)P] 



(4.2) 
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We transform this by defining the action of on a function P' = P/x2 as follows 

= + dl){x 2 P') + -d^xxXzP') + —d 2 [{x% - 1)P'] 

X 2 X 2 X 2 

= di[(d 1 +x 1 )P'] + —d 2 [x 2 (d 2 +x 2 )P'} . (4.3) 

X2 

We now consider T' to be an operator acting in three-dimensional space, with cylindrical 
polar coordinates (r, if, z). We identify r — X2, and z = x-y, and take P' to be a function 
which is restricted so that it has cylindrical symmetry, being independent of ip. With 
this interpretation, we can add differentials with respect to 4> to the definition of J°q, and 
write 

P Q = l -d r [r{d r + r)} + ^3% + d z (d z +z) = V.(x + V) (4.4) 

which is the Fokker-Planck operator for isotropic diffusion in three-dimensional space 
(with D = 1), with a drift velocity v = — x. Thus we have transformed the two- 
dimensional Fokker-Planck equation to a three-dimensional one with a very simple unper- 
turbed velocity field. It is convenient to work with Cartesian coordinates x = (x, y, z) in 
the three-dimensional space, having the usual relation to the cylindrical polar coordinates 
(r, if, z). The Fokker-Planck equation is then 

d v P> = V 2 P' + V.[(r - ev'JP'} = T'P' = [£ + eT x ]P' (4.5) 

where, in Cartesian coordinates, has components 

v' u = —2xz , 

v' 12 = -2yz, (4.6) 
v'i 3 = - 



-z 2 + T(x 2 +y 2 ). 



We now transform the Fokker-Planck T' operator so that is transformed into a very 
simple Hermitian operator, by writing 

H = cxp($ /2)J £ "exp(-$o/2) with $ = \{x 2 + y 2 + z 2 ) . (4.7) 

We find (on writing (x, y, z) — (z\, z 2 , Zs)) 

3 

Ho = exp($ /2)^ exp(-$ /2) - Y.K ~ \ z ) + ^ ( 48 ) 

i=i 

so that H.q is (apart from a negative multiplicative factor) the Hamiltonian operator for 
a three-dimensional harmonic oscillator. The spectrum of Hq is the set of non-positive 
integers (0 ,-1, —2, —3. ..). The eigenfunctions of Ho are generated by raising and lowering 
operators l|Dirac 19 30): 

hi = \zi + d Zi , a+ = \zi - d Zi . (4.9) 
The transformed perturbation operator is 

= 2a+ Zl z 3 + 2a+z 2 z 3 - aj[zl - T(z 2 + z|)] . (4.10) 
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Instead of solving the Fokker-Planck equation jF' P' = we attempt to solve HQ = 0, 
where Q = exp($ /2)P'. 

Now consider how to obtain the Lyapunov exponent from the function Q. We have 
Ai = 76(2:3). We calculate the average of Z3 as follows 







Z3) = / dr dz P(r, z) 



OO f-1-K 

r dr dip / dz zP'(r,z) 



Jo 

00 poo p 00 



/oo poo p 00 
/ / dx dy dz z exp(-<Z> Q /2)Q(x, y, z) . (4.11) 
-oc J —00 J —00 

Now we change the notation, using a variant of the Dirac notation to represent the 
function Q by a 'ket vector' \Q). Allowing for the possibility that \Q) is not normalised, 
we write 

^ = (0000 1 h I Q) _ (0ooo|a3|Q) , 4 
(0ooo|<9) (0000 |Q) 

where |0ooo) is the ground-state eigenfunction of Ho, given by the function exp(— $o/2) = 
exp[-(z2 + z| + z|)/4]. 

We calculate \Q) by perturbation theory: writing 

10) = \Qo) + e\Qi) + e 2 \Q 2 ) + - (4.13) 

we find that the functions \Qk) satisfy the recursion relation 

\Qk+x) = -Ha ^ilQfc) . (4.14) 

At first sight this appears to be ill-defined because one of the eigenvalues of Ho is zero, so 
that the inverse of Ho is only defined for the subspace of functions which are orthogonal 
to the ground state, |0ooo)- However, because all of the components of Hi have a creation 
operator as a left factor, the function Hiltp) is orthogonal to for any function \ip), so 
that (|4.14|l is in fact well-defined. The iteration starts with \Qo) — |0ooo)- The functions 
\Qk) should all have rotational symmetry about the z-axis. The angular- momentum 
operator J$ = p\Z2 — P2Z1 commutes with both Ho and H\ (that is, [Ho, ^3] = and 
[Hi, J3] — where we use square brackets for the commutator, [.4, B] = AB — BA). The 
operators Ho and H\ can therefore be simultaneously reduced to block diagonal form, 
with blocks labelled by eigenvalues of j% . We can restrict ourselves to the subspace where 
the eigenvalue of J-& is zero. The functions of this subspace are generated from the ground 
state using transformed raising and lowering operators, defined as follows: 

Q + = -^=(ai + ia 2 ) , &_ = -^=(61 — i&a) . (4-15) 

The transformed operators a + and a_ satisfy [a±,aj_] = I (where I is the identity 
operator), which is the fundamental relation describing harmonic-oscillator raising and 
lowering operators. Expressing Ho and J% in terms of these operators, we find 

Ji = d_Q!_ — ci^ci + , Ho = — (c^_a- + a\_a + + 0^0,3) . (4-16) 

Using results from lDirac (1930|) . we see that both Ho and S3 are linear combinations of 
harmonic-oscillator Hamiltonians, aLa_, a^a+ and 0,30,3. The nth eigenfunction \4>n) of 
a harmonic-oscillator Hamiltonian a^a is obtained from its ground state |0o) by repeated 
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application of the raising operator a/: 

\cf> n ) = 4^)"!^) (4.17) 



and this eigenfunction has eigenvalue n. Thus eigenfunctions of 7io and ^3 with zero 
angular momentum are constructed as follows 

\^ nm ) = l-L=(aL)"(4)"(4) m |^ooo) (4.18) 

for n = 0, 1, . . . and m = 0, 1, . . .. The corresponding eigenvalues of H.q are — 2n — m. The 
functions \Qk) are expanded in terms of the \ipnm), with coefficients a/nm- 

00 00 

IQ*) = EE 4^1^™) • (4.19) 

n— m— 

By projecting equation (|4.14(l onto the vector |Vw) and using the fact that the eigen- 
vectors |VVm') of H form a complete basis, the iteration can be expressed as follows 
[for (n,m) ^ (0,0)]: 

n'— m'— 

The matrix elements (Vw|7ii|VVm') are readily computed using the algebraic properties 
of the raising and lowering operators, as discussed in lDirac (T930|) . The coefficients aim 
are then calculated recursively. This allows us to obtain the functions \Qk)- The lowest 
order is |Qo) = l^ooo)- Its contribution to Ai vanishes in view of 1)4. 10[) . The leading order 
is 

IQi) = -^11) - l^oi) - ^03) + 2r|^oi) + y |Vn) • (4.21) 

The next order, IQ2), does not contribute to Ai since H.i\Qi) does not contain IV'oi)- In 
fact, only odd orders contain ji/'oi) and thus give non-zero contributions to Ai. We also 
find that the denominator in (|4.13|l is unity at all orders. The final result is: 

00 

Ai = 7 e^ Q (r)e 2 '- 1 (4.22) 

where the first five coefficients ci(T) are 

ci(r) = -l + 2r (4.23) 

c 2 (r) = -5 + 2or- 16 r 2 

c 3 (r) = -60 + 360 r - 568 T 2 + 272 T 3 

c 4 (r) = -1105 + 8840 r - 61936/3 T 2 + 58432/3 T 3 - 19648/3 V 4 
c 5 (r) = -27120 + 271200 T- 7507040/9 T 2 + 3492160/3 T 3 
-2316032/3 T 4 + 1785856/9 T 5 



Th e coefficients in l|4.23|l have a growth which is typical of asymptotic series, as discussed 
bv lDingle (T973|) . Figure ^jp shows approximations to the Lyapunov exponent Ai for 
T = 0.45. Shown are six different partial sums of the series 1)4. 23[) . including terms up 
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0.0 0.5 e 1.0 0.0 0.1 0.2 f 0.3 0.4 0.0 0.2 f 0.4 



Figure 2. a Lyapunov exponent as a function of e: results from eqs. 13,2011 . 113.2111 and 13.221 
are shown as solid lines, those from direct simulations as symbols, T — 1/3 (o), r = 1 (□), 
and r = 2 (o). b Lyapunov exponent \i/(ey) versus e for T — 0.45. Shown are results from 
simulations (o) as well as from the asymptotic series 1)4,230 summed to orders / max = 1, . . . , 6: up 
to e*(7 max ) full lines and from e*(l msx ) dashed lines, c Phase diagram in the e-F-plane. Shown 
are results from numerical simulations of 12. U . o, summation of the asymptotic series l|4.230 
summed to k*(e,T), blue line, as well as results from Langevin simulations (red line). 



to I 

max; with l max — 1, . . . , 6. For a given value of c, there is an optimal value of /max; 
which we term i max , denned by the criterion that the term in 1)4. 23|) with index ^ ax is 
smallest in magnitude. The function Z max (e) can be inverted, its inverse e*(Z max ) is the 
value of e for which the Z max term is optimal. For values of e less than the e*(7 max ) the 
results are shown as solid lines. Beyond this optimal value of e, the results are shown as 
dashed lines. The results show that the series agrees well with the numerical simulation 
up e*, as would be expected for an asymptotic series. 

The phase boundary in the e-F-plane is determined by the condition Ai = 0. Figure^ 
shows this phase boundary as determined from truncating the series l|4.23ll at the optimal 
order Z max (e) (blue line). This asymptotic result is shown for values of e up to « 0.2. 
Beyond this range, the asymptotic approximation becomes increasingly inaccurate. Also 
shown, in the same plot, are results obtained from the Langevin equations [eqs. (|3.20(l . 
H3.21[l and l)3.22[l ]. and from direct simulations (o). The results show that the coalescing 
phase disappears as the effect of inertia is increased: the coalescing disappears altogether 
for the case of pure potential flow (T = |) at e 0.33. 

One notable difference between the t hree-dimensional calculation presented here and 
the two-dimensional case considered in l(Mehlig fc Wilkinson 2 004'! is that in the two- 
dimensional case the phase boundary has a power series in e which vanishes identically 
(and the phase line is therefore non-analytic). 



5. Effect of dispersion of particle masses 

In most naturally occurring aerosols the suspended particles have different mass m, 
and particles of differing sizes will also have different values of the damping coefficient 7. 
It is important to consider whether particles still have a tendency to coalesce even when 
the particles have differing values of m and 7: we argue that the path coalescence effect 
is not destroyed by small mass dispersion. The argument can be adapted to dispersion 
of the damping coefficient, reaching the same conclusion. 

Assume that the path-coalescence effect occurs for particles of mass m. Compare the 
motion of this reference particle with that of an initially nearby particle with mass m+5m. 
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The reference particle has equation of motion 

mx = — 'ymx + f(x, t) . (5-1) 

Writing f(x + 6x,t) = f(x,t) + F(t)5x + 0(5x 2 ), the equation of motion for the other 
particle is 

(to + 5m)(x + 5x) = -7771(1 + 6x) + f(x, t) + F(t)Sx + 0(5x 2 ) . (5.2) 

Collecting the terms which are first order in 5x, we obtain a linearised equation of motion 
for Sx: 

— mSi — jmSx + F{t)5x — 5m x . (5-3) 

This is an inhomogeneous differential equation for the separation 5x between two parti- 
cles, with a driving term proportional to their mass difference 5m. The solution of this 
equation can be constructed from a Green's function satisfying G(t, to) 

-m—- 7m —+F(t)G = S(t-t ) (5.4) 

with G(t,to) = for t < to. The solution of equation (|5.3|l is 

f* d 2 x(t') 
5x{t) = 5m / dt' G{t,t')—-)^ . (5.5) 
Jo dt 

For t > to, equation (|5.4|l is the equation for small displacements of trajectories of parti- 
cles with the same mass. We know that in the path-coalescing phase the solutions have 
a negative value of the Lyapunov exponent Ai, and that they therefore decay exponen- 
tially at large time. In the case where G(t,t') is bounded by an exponentially decreasing 
function, such that |G(t,t')| < v4exp(— Ai|t — 1'|), equation (|5.5|l remains finite as t — * oo. 
For sufficiently large A, the probability of this inequality being violated is extremely 
small. This indicates that in the path-coalescing phase the solution (|5.5|l remains finite 
as t — > oo, except for very rare events. The conclusion is that, when Ai < 0, two initially 
close particles with nearly equal mass will remain in close proximity for a very long time. 



6. Discussion 

In this paper we described the path-coalescence transition, and showed that the tran- 
sition point is determined by the change of sign of a Lyapunov exponent. We showed 
that in general the Lyapunov exponent is determined from an expectation value of a 
variable of a simple dynamical system, equations l|2.13|l . which is driven by stochastic 
forcing functions. We considered the solution of these equations in a particular limiting 
case, where the dimensionless parameters satisfy \ ^ v ^ 1: by mapping the continuous 
differential equations into a pair of coupled Langevin equations. We used these Langevin 
equations to produce a rather complete description of the phase transition in that limit. 
The remainder of these concluding remarks discuss how equations (|2.13|) can be used in 
the case where these inequalities are not satisfied. 

In order to solve these differential equations it is necessary to characterise the statistics 
of the stochastic driving terms F^{t). These terms contain information about the strain- 
rate of the field evaluated at a point along the reference particle trajectory. There are 
two possibilities: 

Case A: the statistics of the strain-rate tensor along a trajectory may be indistinguish- 
able from those sampled along a randomly chosen trajectory. 

Case B: the trajectory may select regions where the strain-rate tensor has atypical 
properties, for example by tracking points where the velocity vector u vanishes. 
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If case A is realised there are two further possibilities: 

Case Al: the trajectory r(t) is sufficiently slowly moving that the displacement over 
time r is small compared to £. In this case the statistics of FL are those of a randomly 
chosen static point, and the correlation time of F^if) will be r. 

Case A2: if the trajectory r(t) is moving sufficiently rapidly that its displacement in 
time t is large compared to £, then the correlation time of F!j(t) will be smaller than 
t because the loss of correlations results primarily from changing the position at which 
dui/dxj(r,t) is sampled. 

The limit which was investigated in detail in this paper (x -C v -C 1) is an example of 
case Al. In cases where \ an d v approach different limits however, all three possibilities 
can occur in the system described by equations (f 1 . 1 1> to Ijl.^JI . 
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